1. Field of the Invention
The present invention relates to computerized simulation of what are known as giant reservoirs.
2. Description of the Related Art
In the oil and gas industries, the development of large underground hydrocarbon reservoirs often requires the building of computer simulation models. These underground hydrocarbon reservoirs are often complex rock formations which contain fluids in the form of petroleum fluid mixtures and water, which exist in two or more fluid phases. During production, the petroleum mixture is produced by wells drilled into and completed in these rock formations. Sometimes, fluids such as water and/or gases are also injected into these rock formations to improve the recovery of the petroleum fluids. Oil and gas companies have come to depend on reservoir simulation as an essential tool to enhance their ability to exploit their petroleum reserves.
Reservoir simulation belongs to the general domain of flow in porous media simulation. However, reservoir simulation normally involves multiple hydrocarbon components and multiple fluid phases in an underground geological formation which is under high pressure and temperature. The chemical phase behavior of these hydrocarbon fluids and the included groundwater has to be taken into account in these simulators.
The simulation models contain the data which describe the specific geometry of the rock formations and the wells, the fluid and rock property data, as well as production and injection history pertaining to the specific reservoirs of the oil or gas field in question. The simulation models are formed by a simulator (known as a reservoir simulator) which is a suite of computer programs run on a data processing system. The reservoir simulator which runs these models is a computer implemented numerical methodology, or coded algorithms and data constructs of an underlying mathematical model. The mathematical model which represents the physics of fluid movements in these hydrocarbon reservoirs is a system of nonlinear partial differential equations which describe the transient multiple-phase, multiple-component fluid flow and material balance behaviors in these reservoirs induced by the production and/or injection of fluids, as well as the pressure-volume-temperature (PVT) relationships of the reservoir fluids.
A reservoir simulator simulates the multiphase multicomponent fluid flow and material balance in subterranean reservoirs and the included surrounding porous rock formations by subdividing the volume into contiguous cells known as grid blocks. A grid block is the basic finite volume where the underlying mathematical model is applied. The number of grid blocks varies depends on the resolution needed for the simulation and the size of the reservoirs in question. For a large reservoir, such as the type known in the industry as a giant reservoir, which may have multi-billion barrels of original oil-in-place (OOIP), the number of grid cells can be in the hundreds of millions to over a billion, in order to have adequate resolution to represent flow dynamics, formation rack porosity and permeability heterogeneity, and many other geologic and depositional complexities within the reservoir. Simulation of this size reservoir can be termed giga-cell reservoir simulation.
Depending on the depositional history of the reservoir layering and the subsequent geological or erosional processes which shaped it, the geometries of the rock layering and hence the description of the grid cells can be very complex. FIG. 1 is a schematic diagram of a typical structural grid G using what is known as corner-point-geometry (CPG) description. The geometric description of grid cells often contains discontinuities such as faults and fractures, pinch-outs and shale barriers. These discontinuities lead to an unstructured cell connectivity requirement even for reservoir simulation models which use a structured grid, such as shown in FIG. 1. These irregularities are known as non-neighbor connections.
Geologic complexity is one driver for reservoir simulation to use unstructured cell connectivities where the flow geometry can be more accurately represented. FIG. 2A is a display of a portion of an example of what is known as an unstructured PEBI grid U for a reservoir of interest with a number of multilateral wells 20 indicated at their respective locations in the reservoir. FIG. 2B is an enlarged view of a portion of FIG. 2A indicated by reference numeral 22. As can be seen in FIG. 2B, the discretization of the grid cells in the vicinity of the multilateral wells 20 in the unstructured grid U honors the local flow geometry rather than being structured and rectilinear as in the case of the structured grid G of FIG. 1.
Another driver is thus the complex well geometry as typical wells drilled today are either horizontal or multilateral wells strategically positioned to maximize reservoir contacts and hydrocarbon recovery. Grids such as those of FIGS. 2A and 2B which honor local flow directions near the wells provide better numerical accuracy. Modern reservoir simulators require solver methodology to be able to address the irregular inter-cell connections. This leads to the requirement to solve large spare linear systems involving general unstructured matrices.
Accurate near-well flow modeling and the modeling of cross-flow is another important aspect in modern reservoir simulation. The solver utilized must include fully-implicit fully-coupled well modeling with the reservoir system in order for the simulation to be accurate and robust. The fully coupled solution must also be fast and scalable in the context of massive parallelism in the computation.
The transient solution of the multiphase multicomponent system involves the evolution of mass and energy conservation in a sequence of time steps from the initial condition of the reservoir. For each time step, the system of nonlinear discrete equations for each finite volume is linearized using what is known as the generalized Newton's method. A general species conservation equation for the component i is given by:
                                                                        ∂                                  c                  i                                                            ∂                t                                      +                          ∇                              ·                                                      ∑                                          j                      =                      1                                                              n                      p                                                        ⁢                                      (                                                                                            ρ                          j                                                ⁢                                                  ω                          ij                                                ⁢                                                  u                          j                                                                    -                                              ϕ                        ⁢                                                                                                  ⁢                                                  ρ                          j                                                ⁢                                                  S                          j                                                ⁢                                                                                                            D                              →                                                        ij                                                    ·                                                      ∇                                                          ω                              ij                                                                                                                                            )                                                                                =                                    ϕ              ⁢                                                ∑                                      j                    =                    1                                                        n                    p                                                  ⁢                                  (                                                            S                      j                                        ⁢                                          R                      ij                                                        )                                                      +                                          (                                  1                  -                  ϕ                                )                            ⁢                              R                is                                      +                                          q                .                            i                                      ⁢                                  ⁢                                  ⁢        where                            (        1        )                                                          ⁢                              c            i                    =                                    ϕ              ⁢                                                ∑                                      j                    =                    1                                                        n                    p                                                  ⁢                                                      ρ                    j                                    ⁢                                      ω                    ij                                    ⁢                                      S                    j                                                                        +                                          (                                  1                  -                  ϕ                                )                            ⁢                              ρ                s                            ⁢                              ω                is                                                                        (        2        )            
If dispersion, chemical reaction and absorption are ignored, the species equation simplifies to:
                                                        1                              V                b                                      ⁢                                          ∂                                  n                  i                  t                                                            ∂                t                                              +                      ∇                          ·                                                ∑                                      j                    =                    1                                                        n                    p                                                  ⁢                                  (                                                            ρ                      j                                        ⁢                                          x                      ij                                        ⁢                                          u                      j                                                        )                                                                    =                              q            .                    i                                    (        3        )            
Since the pore space of porous medium must be filled with fluids present, the pore volume must be equal to the total fluid volume. This can be expressed as:
                                          ∑                          i              =              1                                      n              p                                ⁢                      V            j                          =                  V          ϕ                                    (        4        )            where the pore volume Vφ is a function of pressure alone and described as:
                              V          ϕ                =                              V            ϕ            ref                    ⁢                      ⅇ                                          C                r                            (                                                P                  ref                                -                P                            )                                                          (        5        )            
Pressure and the overall number of moles are the primary variables. For closure, the other equations used are:
                              Constraint          ⁢                                          ⁢          on          ⁢                                          ⁢          mole          ⁢                                          ⁢          fractions          ⁢                                          ⁢          for          ⁢                                          ⁢          each          ⁢                                          ⁢          phase          ⁢                      :                    ⁢                                          ⁢                                    ∑                              i                =                1                                            n                c                                      ⁢                          x              ij                                      =        1                            (        6        )                                                      Constraint            ⁢                                                  ⁢            on            ⁢                                                  ⁢            total            ⁢                                                  ⁢            moles            ⁢                                                  ⁢            per            ⁢                                                  ⁢            component            ⁢                          :                        ⁢                                                  ⁢                                          ∑                                  j                  =                  1                                                  n                  p                                            ⁢                                                n                  j                  p                                ⁢                                  x                  ij                                                              =                      n            i            t                          ⁢                                                      (        7        )                                          Constraint          ⁢                                          ⁢          on          ⁢                                          ⁢          fluid          ⁢                                          ⁢          saturations          ⁢                      :                    ⁢                                          ⁢                                    ∑                              i                =                1                                            n                p                                      ⁢                          S              j                                      =        1                            (        8        )                                          wherein          ⁢                      :                    ⁢                                          ⁢                      S            j                          =                              V            j                                              ∑                              j                =                1                                            n                p                                      ⁢                          V              j                                                          (        9        )                                          and          ⁢                                          ⁢                      V            j                          =                              n            j            p                                ρ            j                                              (        10        )            
Phase velocities are described by Darcy's Law:uj=−Kλj(∇Pj−γj∇D)   (11)where K is the permeability tensor defined as:
                    K        =                  [                                                                      k                  xx                                                                              k                  xy                                                                              k                  xz                                                                                                      k                  yx                                                                              k                  yy                                                                              k                  yz                                                                                                      k                  zx                                                                              k                  zy                                                                              k                  zz                                                              ]                                    (        12        )            
The symbols in the above equations have these meanings:
Symbolppressureqproduction ratexiMole fractionVjPhase VolumeSjPhase SaturationciOverall Concentration of species iφporosityρdensityμviscosityωmass fractionRHomogeneous reaction rateDDispersion CoefficientuvelocityVφRock pore volumenitOverall number of mole
Subscriptsicomponent indexjphase index
Superscripts:refreferencepa fluid phasettotal
In the industry, this is referred to as Newtonian iteration or nonlinear iteration. At each Newtonian iteration, a linear system of equations is constructed where the matrix, known as the Jacobian matrix, and the right-hand-side vector, known as the residuals are used to solve for the changes in the primary variables of the system.
Current industry practice for solving the linear system of equations is via a preconditioned iterative method. The iterative algorithm can be one of three forms: one known as the ORTHOMIN; the second known as the GMRES, or the third known as the BICGSTAB method. For robustness, the preconditioner used, rather than the choice of the iterators, has been more important. Typical state-of-the-art preconditioners used in the reservoir simulation industry can be the one-level preconditioner such as the nested factorization (NF) method, or variants of the incomplete LU factorization (ILU) method. These are widely practiced in today's commercial simulators.
More recent practice involves the use of multi-stage preconditioning. The nature of the governing equations in reservoir simulation is mixed parabolic-hyperbolic. For treating the parabolic part of the system, the multi-level algebraic multigrid (AMG) or geometric multigrid (GMG) can be very effective if the pressure components of the original problem are first decoupled from the full system via what is known as the constraint pressure residual (CPR) algorithm. The approximate pressure problem can be attacked effectively via a multi-level method such as AMG or GMG. The remaining hyperbolic component of the residuals can be readily dealt with using a suitable and cheaper one-level preconditioner. This in effect is what is known as a divide-and-conquer preconditioning strategy, and is known in the art. However, the GMG multi-level method cannot, so far as is known, handle non-neighbor connections. Thus, the AMG method is the generalized multi-level method typically used to accommodate realistic simulation problems.
The measure of efficiency of a computational processing algorithm in parallel computing is its scalability. A method is perfectly scalable, or has a 100 percent parallel efficiency, if it takes one hour to solve a computational problem on the computer with a single CPU, and then if the work is exactly divided out and given to two CPU's, the time to solve the same problem is 0.5 hour, and the time to solve the same problem using four CPU's is 0.25 hour and so on. That is, there is no parallelization overhead. A perfectly scalable method would be an ideal situation. In reality, many reasons can cause the solution time on a real system to be far from this ideal.
The practices or strategies described above work very well for serial computing, using a single core of a single CPU, and for small scale parallelism using tens of CPU cores. However, they suffer significant loss of parallel scalability as the number of processes increases. This is because the nested factorization (NF) and incomplete LU factorization (ILU) described above have significant recursive components which cannot be easily parallelized. In addition, multi-level methods such as algebraic multigrid (AMG) incur significant communication overhead in the aggregation and smoothing at coarsened grid levels which reduces parallel scalability.
Increasingly, simulation models of reservoirs and oil/gas fields have become more large and complex. At the same time, computer hardware has evolved rapidly to become inexpensive and fast. The current trend for hardware evolution is towards massive parallelism: from high performance computing or HPC system with single-core CPU's in the recent past, to currently used PC clusters composed of thousands of multi-core CPU's, to the recent advent of the heterogeneous CPU-GPU system. Simulation practitioners can now use inexpensive hardware with tens of thousands of compute cores working concurrently to solve a large-scale simulation problem.
So far as is known, existing solver methods did not handle giga-scale reservoir simulation due to poor scalability issues and/or robustness issues. Presently available solver methodology is limited to about ten million grid cells and parallel scalability in the tens of processes. Without massively parallel scalability, a simulation run for a history match simulation or a prediction simulation becomes too slow and impractical. This limits the size of the models which can be used in flow simulation. For the large reservoirs, this means excessive upscaling.
Several other factors about the reservoir further increased the complexity of giga-scale reservoir simulation. The treatment of common geologic discontinuities such as fractures, faults, pinch-outs, and shale barriers posed difficulties for robust scalable solver methods. The complex flow geometry near multilateral wells was not accurately modeled using structured grids. Unstructured grids required specialized unstructured methods and made the computational load harder to balance among processors and harder to solve robustly and efficiently.
Also, the fully-coupled fully-implicit well solution is a required method for accurate reservoir simulation. This is another factor which can degrade scalability, especially when the reservoir contains thousands of complex multilateral wells. Further, many petroleum reservoir rocks have multi-modal porosity-permeability and multi-scale fracture networks.